import numpy as np
import csv
import random
from datetime import datetime
from time import strftime
import itertools as itertools
import math as math

'''
revision of test25 for Mar21
to calculate the false nagetive, false positive

'''

def getCurTime():
    """
    get current time
    Return value of the date string format(%Y-%m-%d %H:%M:%S)
    """
    format='%Y-%m-%d %H:%M:%S'
    sdate = None
    cdate = datetime.now()
    try:
        sdate = cdate.strftime(format)
    except:
        raise ValueError
    return sdate

def build_unit_list(inputCSV):
    sKey = np.array([])
    fn = inputCSV
    ra = csv.DictReader(file(fn), dialect="excel")
    i = 0
    for record in ra:
        #print record[ra.fieldnames[0]], type(record[ra.fieldnames[-1]])
        temp = [int(float(record["GIST_ID"])), int(float(record[cancerFN])),int(float(record[popFN])), 0]
        sKey =np.append(sKey, temp)
        i += 1
    sKey.shape=(-1,4)
    return sKey

def build_region_list(inputCSV):
    fn = inputCSV
    ra = csv.DictReader(file(fn), dialect="excel")
    temp_list = np.array([])
    i = 0
    for record in ra:
        #print record[ra.fieldnames[0]], type(record[ra.fieldnames[-1]])
        temp_list = np.append(temp_list, [int(float(record[ra.fieldnames[0]])), int(float(record[ra.fieldnames[-1]]))])
        #unit_attri[i,3] = int(float(record[ra.fieldnames[-1]]))
        i = i + 1
    temp_list.shape = (-1,2)
    return temp_list
        
def cal_bernoulli_likelihood_ratio(iObZ, iObT, iPopZ, iPopT):
    # iObZ: the observed number in the certain Zone, e.g. the number of cancer cases in a region
    # iObT: the total observed number, e.g. the total number of cancer cases
    # iPopZ: the population in the certain Zone 
    # iPopT: the total population
    # return the natural logarithm of ratio
    temp_lambda = 0.0
    temp_z_p = (0.0 + iObZ)/iPopZ
    temp_out_p = (0.0 + iObT - iObZ)/(iPopT - iPopZ)
    temp_t_p = (0.0 + iObT)/iPopT
    if temp_z_p > temp_out_p:
        #temp_denominator = temp_t_p ** int(iObT)
        #temp_numerator = (temp_z_p**int(iObZ))*(temp_out_p**int(iObT - iObZ))
        #temp_lambda = temp_numerator/temp_denominator
        temp_lambda = iObZ * math.log(temp_z_p) + (iPopZ - iObZ) * math.log(1 - temp_z_p) +\
                        (iObT - iObZ)* math.log(temp_out_p) + ((iPopT - iPopZ) - (iObT - iObZ))* math.log(1 - temp_out_p) \
                        - iObT * math.log(temp_t_p) - (iPopT - iObT) * math.log(1 - temp_t_p)
    return temp_lambda

def cal_poisson_likelihood_ratio(iObZ, iObT, iPopZ, iPopT):
    # iObZ: the observed number in the certain Zone, e.g. the number of cancer cases in a region
    # iObT: the total observed number, e.g. the total number of cancer cases
    # iPopZ: the population in the certain Zone 
    # iPopT: the total population
    # return the natural logarithm of ratio
    temp_lambda = 0.0
    temp_z_p = (0.0 + iObZ)/iPopZ
    temp_out_p = (0.0 + iObT - iObZ)/(iPopT - iPopZ)
    temp_t_p = (0.0 + iObT)/iPopT
    if temp_z_p > temp_out_p:
        #temp_denominator = temp_t_p ** int(iObT)
        #temp_numerator = (temp_z_p**int(iObZ))*(temp_out_p**int(iObT - iObZ))
        #temp_lambda = temp_numerator/temp_denominator
        temp_lambda = iObZ * math.log(temp_z_p) + (iObT - iObZ)* math.log(temp_out_p) - iObT * math.log(temp_t_p)
    return temp_lambda

def cal_region_likelihood_ratio(temp_region_attri, temp_total_caner, temp_total_pop):
    # t_unit_attri: the attribution (id, cancer, pop, region_id) of the whole study area
    # t_region_attri: region id for each unit
    
    temp_poisson_ratio = np.array([])
    temp_bernoulli_ratio = np.array([])
    
    for item in temp_region_attri:
        temp_bernoulli_ratio = np.append(temp_bernoulli_ratio, cal_bernoulli_likelihood_ratio(item[0], temp_total_caner, item[1], temp_total_pop))
        temp_poisson_ratio = np.append(temp_poisson_ratio, cal_poisson_likelihood_ratio(item[0], temp_total_caner, item[1], temp_total_pop))
    
    return temp_poisson_ratio, temp_bernoulli_ratio
    #return temp_poisson_ratio

def cal_region_attri(temp_unit_attri, temp_region_id):
    temp_region_attri = np.zeros((len(temp_region_id),3))    #[cancer, pop, count]
    for item in temp_unit_attri:
        temp_region_attri[int(item[3]), 0] += item[1]
        temp_region_attri[int(item[3]), 1] += item[2]
        temp_region_attri[int(item[3]), 2] += 1
    return temp_region_attri
    
def cal_region_cancer(temp_unit_attri, temp_region_id):
    temp_region_cancer = np.zeros(len(temp_region_id))
    for item in temp_unit_attri:
        temp_region_cancer[int(item[3])] += item[1]
    return temp_region_cancer

def cal_region_cancer_d(temp_cancer_threshold, temp_total_caner, temp_total_pop):
    temp_cancer = np.zeros(len(temp_cancer_threshold))
    i = 0
    while i < temp_total_caner:
        temp = np.where(temp_cancer_threshold < random.randint(0, temp_total_pop-1))
        temp_cancer[len(temp[0])] += 1
        '''
        temp = random.randint(0, temp_total_pop-1)
        temp_left = 0
        for k in range(len(temp_cancer_threshold)):
            if (temp >= temp_left) & (temp_cancer_threshold[k] > temp):
                temp_id = k
                break
            temp_left = temp_cancer_threshold[k]
        temp_cancer[k] += 1
        '''
        i += 1
    return temp_cancer

#--------------------------------------------------------------------------
#MAIN

if __name__ == "__main__":
    print "begin at " + getCurTime()

    filePath = "C:/_DATA/CancerData/"
    unitCSV = "C:/_DATA/CancerData/LATrtCancerS.csv"
    popFN = "POP_2000"
    shuffle_id = 1  
    # 0: shuffle the cancer count for each unit, 1: random generate the cancer count for each unit conditioned on the total count
    
    for year in range(2,5):
        print "calculating 200" + str(year) + "..."

        if year == 2:
            soMethod = "CLK"
        elif year == 3:
            soMethod = "ALK"
        elif year == 4:
            soMethod = "WARD"
        else:
            print "error in soMethod!"
            raise SystemExit
        regionCSV = "C:/_DATA/CancerData/test/Apr1_real_data/total5_"+soMethod+"_SO_POP_min_20000.csv"
        
        # COUNT02    COUNT03    COUNT04    COUNT05    COUNT06    COUNT02_06
        cancerFN = "COUNT02_06"

        
        unit_attri = build_unit_list(unitCSV)  #unit_attri: [id, cancer, pop, region_id]
        total_caner = np.sum(unit_attri[:,1])
        total_pop = np.sum(unit_attri[:,2])
        region_id = build_region_list(regionCSV)
    
        for unit in unit_attri:
            for region in region_id:
                if int(unit[0]) == int(region[0]):
                    unit[3] = region[1]
                    break
        region_id = np.unique(unit_attri[:,3])
        region_id = np.sort(region_id)
        
        region_attri = cal_region_attri(unit_attri, region_id)  #region_attri: [cancer, pop, count]
        
        poisson_likelihood_ratio, bernoulli_likelihood_ratio = cal_region_likelihood_ratio(region_attri, total_caner, total_pop)
        #poisson_likelihood_ratio = cal_region_likelihood_ratio(region_attri, total_caner, total_pop)
        if shuffle_id == 1:
            cancer_threshold = np.array([])
            temp = 0
            for item in region_attri:
                cancer_threshold = np.append(cancer_threshold, temp + item[1])
                temp = cancer_threshold[-1]
    
        '''
        if shuffle_id == 1:
            cancer_data_a = np.array([])
            cancer_data_b = np.array([])
            for item in region_attri:
                cancer_data = np.append(cancer_data, np.ones(int(item[0])))
                #cancer_data = np.append(cancer_data, np.zeros(int(item[1])))
        '''
        print "Permutation started at " +  getCurTime() 
    
        repeat = 9999
        iCnt = 0
        iTenth = 1
        fInterval = 2000
        for i in range(repeat):
            if shuffle_id == 0:
                random.shuffle(unit_attri[:,1])
                region_attri[:,0] = cal_region_cancer(unit_attri, region_id)
            elif shuffle_id == 1:
                #random.shuffle(cancer_data)
                region_attri[:,0] = cal_region_cancer_d(cancer_threshold, total_caner, total_pop)
            else:
                print "error!"
                raise SystemExit
            temp_poisson, temp_bernoulli = cal_region_likelihood_ratio(region_attri, total_caner, total_pop)
            #temp_poisson = cal_region_likelihood_ratio(region_attri, total_caner, total_pop)
            poisson_likelihood_ratio = np.append(poisson_likelihood_ratio, temp_poisson)
            bernoulli_likelihood_ratio = np.append(bernoulli_likelihood_ratio, temp_bernoulli)
            iCnt += 1
            if iCnt > fInterval:
                print str(iTenth*2) + "0% permutation finished at " + getCurTime()
                iTenth += 1
                iCnt = 0
        # print likelihood_ratio
        poisson_likelihood_ratio.shape = (repeat+1, -1)
        bernoulli_likelihood_ratio.shape = (repeat+1, -1)
        #np.savetxt("C:/temp/likelihood_ratio.csv", likelihood_ratio, delimiter=',')
        
        pvalue = np.zeros((len(region_id),3))
        iLen = poisson_likelihood_ratio.shape
        i = 0
        while i < iLen[1]:
            pvalue[i, 1] = 1.0
            if poisson_likelihood_ratio[0,i] > 0:
                temp = np.where(poisson_likelihood_ratio[:,i] > poisson_likelihood_ratio[0,i])
                pvalue[i, 1] = len(temp[0])/(repeat+1.0)
            pvalue[i, 0] = region_id[i]        
            i += 1
            
        
        i = 0
        while i < iLen[1]:
            pvalue[i, 2] = 1.0
            if bernoulli_likelihood_ratio[0,i] > 0:
                temp = np.where(bernoulli_likelihood_ratio[:,i] > bernoulli_likelihood_ratio[0,i])
                pvalue[i, 2] = len(temp[0])/(repeat+1.0)
            i += 1
        
        pvalue_path = regionCSV[:-21] + "_pvalue_1.csv"
        np.savetxt(pvalue_path, pvalue, delimiter=',')
    print "end at " + getCurTime()
    print "========================================================================"  

           
